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Abstract 

We deal with the problem of the final sign of the neutrino asymmetry gen- 
erated by active-sterile neutrino oscillations in the Early Universe solving the 
full momentum dependent quantum kinetic equations. We study the param- 
eter region 10~^ ~ \5m'^\/eV'^ < 10^. For a large range of sin^ 2^0 values the 
sign of the neutrino asymmetry is fixed and does not oscillate. For values of 
mixing parameters in the region 10^^ ~ sin^ 29q -^^ 3 x 10~^ (eV^/|5m^|), the 
neutrino asymmetry appears to undergo rapid oscillations during the period 
where the exponential growth occurs. Our numerical results indicate that the 
oscillations are able to change the neutrino asymmetry sign. The sensitivity 
of the solutions and in particular of the final sign of lepton number to small 
changes in the initial conditions depends whether the number of oscillations 
is high enough. It is however not possible to conclude whether this effect is 
induced by the presence of a numerical error or is an intrinsic feature. As 
the amplitude of the statistical fluctuations is much lower than the numerical 
error, our numerical analysis cannot demonstrate the possibility of a chaoti- 
cal generation of lepton domains. In any case this possibility is confined to a 
special region in the space of mixing parameters and it cannot spoil the com- 
patibility of the ^ Vg solution to the neutrino atmospheric data obtained 
assuming a small mixing of the with an eV — r neutrino. 
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I. INTRODUCTION 



If light sterile neutrinos exist, then this will lead to important implications for early 
Universe cosmology. This is because ordinary-sterile neutrino oscillations generate large 
neutrino asymmetries for the large range of parameters, ~ 10^^ eV'^, sin^ 26o ~ 10^^° 
This is a generic feature of ordinary-sterile neutrino oscillations for this parameter 
range. For ~ 10~^ eV"^, the evolution of neutrino asymmetry is qualitatively quite 
different as collisions are so infrequent and a large neutrino asymmetry cannot be generated 
Interestingly, some people do not currently accept that large neutrino asymmetry is 
generated in the early Universe |10[- We will comment briefly on this later in the paper. 

An important issue which has yet to be fully addressed is the sign of this asymmetry. Is it 
always fixed or can it be random? This is an important issue because a random asymmetry 
may lead to domains with lepton number of different signs . If such domains exist then this 
may lead to observable consequences. For example, active neutrinos crossing the boundaries 
of these lepton domains could undergo a MSW resonance which would lead to a new avenue 
of sterile neutrino production |n]. In Refs. P"^|-P^ the issue of sign of the asymmetry was 
discussed in the approximation that all of the neutrinos have the same momentum (i.e. 
p = 3.15T instead of the Fermi-Dirac distribution). This approximation is not suitable 
for discussing the temperature region where the exponential growth in neutrino asymmetry 
occurs. The reason is that in the average momentum toy-model, all of the neutrinos enter 
the MSW resonance at the same time which significantly enhances the rate at which neutrino 
asymmetry is created at T = T^. The rapid creation of neutrino asymmetry significantly 
reduces the region where the oscillations are adiabatic 0. 

Thus it is clear that the neutrino momentum dependence must be properly taken into 
account. This was done in Ref. where an approximate solution to the quantum kinetic 
equations was derived. This approximate solution was called the 'static approximation' and 
was re-derived in a different way in Ref. ||^ where it was shown that this approximation was 
just the adiabatic limit of the quantum kinetic equations (QKE's) in the region where lepton 
number generation is dominated by collisions. Anyway, in the limit where this approximation 
is valid, it was shown in Ref. that the sign is completely fixed. The static approximation 
is valid for a large range of parameters but is not valid for large sin^ 26o ~ 10~^, 6m? ~ 
—10 eV"^. It breaks down in this region because the neutrino asymmetry is generated so 
rapidly during the exponential growth phase that the quantum kinetic equations are no 
longer adiabatic. Thus, while Ref. partially answers the question of sign, it does not give 
the complete answer. The purpose of this paper is to examine the issue of the sign of the 
asymmetry by numerically solving the quantum kinetic equations. 

The outline of this paper is as follows: In section 2 we present some necessary prelim- 
inary discussion on active-sterile neutrino oscillations in the early Universe. In section 3 
we examine the likely size of the statistical fluctuations in the early Universe. In section 
4 we describe the numerical results of our study of the region of parameter space where 
the sign of the neutrino asymmetry is fixed. Using the results of section 3, we are able 
to conclude that in this region the statistical fluctuations cannot have any effect and the 
generated lepton number would have the same sign in all the points of space. In section 5 we 
describe the features of the transition from the region with no oscillations to one where the 
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neutrino asymmetry oscillates for a short period during the exponential growth. In section 
6 we conclude. Also included is an appendix giving some numerical details, which we hope 
will be useful to other workers in the field such as the authors of Ref. IlIO . 



II. PRELIMINARY DISCUSSION 

Our notation/convention for ordinary-sterile neutrino two state mixing is as follows. The 
weak eigenstates Ua {a = e, fi or r) and Ug are linear combinations of two mass eigenstates 
Ua and z/fe, 

Ua = cos OoUa + siu OoUb, Ug = - siu OoUa + COS (1) 

where is the vacuum mixing angle. We define so that cos 2^0 > and we adopt the 
convention that = ml — ml. Recall that the a-type neutrino asymmetry is defined by 

In the above equation, is the number density of photons. Note that when we refer to 
"neutrinos", sometimes we will mean neutrinos and/or antineutrinos. We hope the correct 
meaning will be clear from context. Also, if neutrinos are Majorana particles, then techni- 
cally they are their own antiparticle. Thus, when we refer to "antineutrinos" we obviously 
mean the right-handed helicity state in this case. 

The density matrix for an ordinary neutrino, Ua = e,/i, r), of momentum p 

oscillating with a sterile neutrino in the early Universe can be parameterized as follows: 

Pa/3'(P) = + P(P) • Pa/?'(P) = + P(P) ■ 0^], (3) 

where / is the 2x2 identity matrix, the "polarisation vector" P(p) = Px{p)'i^+Py{p)y+Pz{p)i- 
and a = o"a;X + ayy + cx^z, with cxj being the Pauli matrices. It will be understood that the 
density matrices and the quantities Pi{p) also depend on time t or, equivalently, temperature 
T. The time-temperature relation for me T ~ m^ is dt/dT ~ — Mp/5.44T^, where 
Mp ~ 1.22 X 10^^ MeV is the Planck mass. 

We will normalise the density matrices so that the momentum distributions of Uaip) and 
UsIp) are given by 

fuAp) = \[Poip) + PMUoip), fuM = liPoip) - PMUoip), (4) 

where 

Up) ^ (5) 

1 + exp j 

is the Fermi-Dirac distribution (with zero chemical potential). Similar expressions pertain 
to antineutrinos (with P{p) P{p) and Pq Pq). The evolution of P{p),Po{p) (or 
P{p),Pq{p)) are governed by the equations [p!5|-[r7|j5[1 
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dP dP 

— = V(x) X P(a;) - D{x)[P^{x)S^ + Py{x)y] + ^ z, 



dt 



r(x) 



(6) 



where D{x) = r{x)/2 and T[x) is the total colhsion rate of the weak eigenstate neutrino of 
adimensional momentum x = p/T with the background plasmaQ and feq{x) is the Fermi- 
Dirac distribution: 



1 



1 + exp(x — /ic 



(7) 



where fia= fJ'Ua/T- For anti-neutrinos, i2'a^fJ'a= The chemical potentials /ij,^, /Xp^, 

depend on the neutrino asymmetry. In general, for a distribution in thermal equilibrium 



x'^dx 



4C(3)^ 1 + e 



X-fla 



'tic dj^jc 



4C(3) JO 1 -|- e^'-M5 



(8) 



where C(3) ^ 1.202 is the Riemann zeta function of 3. Expanding out the above equation. 



24C(3) 



7r^(/^a - ^^a) + 6(/i„ - ^^^) ln2 + - /i, 



(9) 



which is an exact equation for /^a= — /^a, otherwise it holds to a good approximation 
provided that 1- For T ~ T^^^ (where T|,^ ^ 2.5 MeV and T^^jJ ^ 3.5 MeV are the 

chemical decoupling temperatures), /i^^ ~ ~/^Pa because processes such as Ua + T^a ^ + e~ 



are rapid enough to make /^q, + /^^ ~ /^e+ + /^e-— 0. However, for IMeV ^ T ^ T, 



dec^ 



weak interactions are rapid enough to approximately thermalise the neutrino momentum 
distributions, but not rapid enough to keep the neutrinos in chemical equilibrium [|. In this 
case, the value of /^q, is approximately frozen at T ~ TJ^^ (taking for definiteness L^^ > 0), 
while the (negative) anti-neutrino chemical potential fJ^a continues decreasing until T ~ 1 
MeV. 

The quantity V(x) is given by |[TB|JT^ 



V(x) = /3(x)x + A(x)z, 



where P{x) and X{x) are 

2xT 



sin 26*0, A(x) 



2^T 



[cos 26*0 — b{x) ± a{x)], 



(10) 



(11) 



^ From Ref. |T^,5[ it is given by r(x) = yG'^pT^x where y ~ 1.27 for v = and y ~ 0.92 for 
= v^, Vr (for me ~ r ~ m^). 

^The chemical and thermal decoupling temperatures are so different because the inelastic collision 
rates are much less than the elastic collision rates. See e.g. Ref. [|l^] for a list of the collision rates. 
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in which the +(— ) sign corresponds to neutrino (anti- neutrino) oscillations. The dimen- 
sionless variables a{x) and b{x) contain the matter effects JTOjl (more precisely they are the 
matter potential divided by 5m?/2xT). For z/q, — > Ug oscillations a(x),6(x) are given by 

10,11 



-4C(3)V2GFTU(")a: _ -4C(3) V2G^rM„x- 



where Gi? is the Fermi constant, My/ is the ly— boson mass, A^. ^ 17 and A^^^ ^ 4.9 (for 
777,g ^ T ~ m^). The important quantity L^") is given by the following expression: 

L(") = + L,^ + L,^ + L,, + = 2L,„ + L, (13) 

where is a small term due to the asymmetry of the electrons and nucleons and is expected 
to be very small, |?7| ~ 5 x 10~^°. We will refer to L*^°^ as the effective total lepton number 
(for the a-neutrino species) since it really needs a name. Note that the quantity L is 
independent of L^^ and its value is currently unknown, but could presumably be calculated 
within a model of baryo-leptogenesisQ. For our numerical work we took two different values 
L = 5 X 10~^° and L = 5 x 10~^^; in this way we could verify that the final value of the total 
lepton number does not depend on this particular choice. The neutrino asymmetry L,^^, on 
the other hand, evolves dynamically and is a function of P, Pq and P, Pq. The equations 
Eq.(^ constitute a closed set of differential equations for the 8 distributions P,Po,P,-Po- 
The 'initial' conditions (for T oo) are simply P^, = Pj^ = 0, Pq = P^ = 1 , assuming here 
for simplicity that initially L,,^ = (and similarly for antiparticles). Note that the vanishing 
of Px,y is just an example of the quantum Zeno effect. 

From a computational point of view it is useful to write down an explicit differential 
equation for the neutrino asymmetry to be solved together (see also [^]). This can be done 
considering that dL^^/dt = —dLj,Jdt and, using the Eq.(^), one easily gets: 



dLi, 

Uct 



dt 8C(3) 



j l3{x)[Py{x) - Pyix)]fo{x) x^dx. (14) 



We mainly employ the useful time saving approximation of integrating the oscillation equa- 
tions, Eq.(|^) [and obviously also Eq.(|l^)] in the region around the MSW resonances. A 
detailed description of the numerical procedures is given in the Appendix. Actually away 
from the resonance the oscillations are typically suppressed by the matter effects or by 
sin^26'o (or both). Thus, this should be a good approximation, which we carefully checked 
by taking larger slices of momentum space around the resonance (an example of this check 
is given in the Appendix). 

Before we finish this section we would like to briefly comment on the static approximation 
The static approximation can be derived by solving the Quantum kinetic equations in 
the adiabatic limit and assuming that the evolution is dominated by collisions. The resulting 



^Alternatively, it maybe set by divine intervention. 
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equations are equivalent to keeping the differential equations for Pq,Pz,Po, Pz [contained in 
Eq.(^] which are functions of Py, Py which are given by 

and Py is similar to the above equation except for the obvious replacement of Pz Pz 
and a —a in the definition of A. The evolution of Ly^ is obtained from Eq.(|l4D above. 
Numerically, the static approximation is much easier to solve, but it is not always a valid 
approximation. It is not valid at low temperatures where coherent MSW transitions are 
important. Also, if the neutrino asymmetry is created fast enough during the exponential 
growth phase then it may not always be valid (which is the case for ~ — 10 eV'^ and 
sin^ 2^0 ~ 10~^, for example). In figures 1-3, we give some numerical examples. For fig. 
la we take Vr ^ oscillations for the parameter choices = —10 eV"^, sin^ 2^o = 10^''. 
Shown is the evolution of the effective total lepton number, L^'^\ with the solid line rep- 
resenting the numerical solution of the quantum kinetic equations integrating around the 
resonances The dashed line is the static approximation discussed above integrating also 
around the resonances. Observe that there is exact agreement between the QKEs and the 
static approximation when both are integrated around the resonances, except at low tem- 
peratures where the static approximation is not valid, as the evolution there is dominated by 
coherent oscillations (MSW effect) as discussed in detail in Ref. 0. In fig. lb we compare 
the static approximation integrating around the resonance with the static approximation 
integrating over the entire momentum range (which we approximate to 0.01 < p/T < 20). 
There is some differences between the static approximation integrated over the entire mo- 
mentum range with the static approximation (or QKE's) integrated around the resonances. 
The difference at high temperature can be qualitatively understood because the resonances 
at high temperature occur at very low momentum, Pres/T ~ 1 where there are very few 
neutrinos. Thus, the region around Pres/T ~ 2 can be somewhat important despite being 
away from the resonance because of the larger number of neutrinos. Also note that the point 
where the exponential growth occurs is modified slightly, but the rate of growth (the slope) 
is approximately unmodified. Whether or not there are oscillations of sign depends on the 
rate at which lepton number is being created during the exponential growth (small changes 
in the value of would not matter). For this reason, the approximation of integrating 
around the resonance should be an acceptable approximation. 

Figures 2a,b are the same as figure la,b except for the different choice of parameters, 
= -100 eV^, sin2 2^o = 10"^. 

As a final check of that static approximation we numerically solve the QKEs for the 
entire momentum range (0.1 < p/T < 12.0) for one example. This is quite CPU time 
consuming and it is most economically done for small 5m^. We take = —0.01 eV'^ 
and sin^ 29o = 10^^. For these parameters, figure 3 compares the numerical solution of the 
QKE with the static approximation, both integrated over the entire momentum range. Note 



^Note that at low temperatures T ~ 5 MeV the repopulation must be taken into account over 
the entire momentum range, 0.01 < p/T < 20. 
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that only the high temperature evolution is shown. As the figure shows we obtain excellent 
agreement as we expect f\ 

Interestingly, Ref. |T^ appears to have (almost) rederived the static approximation (they 
neglect the chemical potentials in the repopulation, which is important for the issue of sign). 
They seem to get a quite different numerical solution when they numerically solve the static 
approximation, and are not able to obtain any solution at all when they numerically solve 
the QKE's. In our opinion, the remarkable agreement between our solution of the static 
approximation and the QKE's (as figures 1-3 illustrate) is a convincing check that we have 
done the numerical work competently. Of course this check has been done previously, and 
there are examples already given in the literature [Q, but since this work was ignored in 
Ref. [llOl, we have taken the trouble to emphasise it again here. 



III. STATISTICAL FLUCTUATIONS 

In this section we want to provide a simple estimation of the statistical fiuctuations in the 
effective total lepton number. These arise simply because of the fluctuations in the number 
of particles within the region, around each oscillating neutrino wavepacket, that determines 
the properties of the medium through the effective potential. 

Let us flrst proceed without specifying the size of this region but we write it in units of 
the interaction length £int = T"^ of the oscillating neutrino with x ~ 2.2 (the peak of the 
distribution). The number of photons within a region of size i at the temperature T is given 
by: 

Since the number of neutrinos is of the same order as the number of photons, it follows that 
the statistical fluctuation of lepton number is of order: 

3 
' 2 



AL = ^ ~ 10-^^ ^ _ . (17) 

Note that since ub ~ 10~^°n^, it follows that A?7 = ^Jn^jn-^ ~ lO^^AL and thus AL^"^ ~ 
AL. 

We are interested to evaluate this fluctuation at the critical temperature T^, when the 
lepton number starts to be generated. The critical temperature can be put in the following 
approximated form 

1 

18 MeV. (18) 




^In the appendix we give further details about the approximation of integrating around the reso- 
nance. We also compare the two different numerical procedures employed. 
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We can thus more usefully express the temperature in units of Tc in Eq. ([T7|), obtaining: 




(19) 



Now what should we take for the size We will present a heuristic argument that 
^ 10~'^iint- First of all, we find numerically that neutrino asymmetries are typically not 



significantly created on time scales less than of order 10~^£j„( i.e. ^int^ ~ C(IO^L) for 
the parameter space of interest (this is true even in the exponential growth phase). This 
means that we only need to consider time scales greater than 10~^£j„t. Furthermore since 
neutrinos are free streaming on such small time scales it follows that the distance scale for the 
creation of L is greater than 10~^£j„t. Thus the characteristic volume relevant for statistical 
fluctuations is larger than or of order (10~^£inf)^. Hence, for \5m?\ < 1000 eV^ we conclude 
that: 

AL(°) ~ O (iQ-^^) . (20) 



In the next sections we will have to find out whether such fluctuations are able to induce a 
random sign of the final lepton number. Note that in addition to statistical fluctuations, it is 
also possible to envisage inhomogeneities which may arise from a model of lepto-baryogenesis 
for example. Such inhomogeneities may potentially be much larger than the statistical fluctu- 
ations but are of course model dependent. Note that the possibility of such inhomogeneities 
has recently been considered in Ref. where it is shown how the diffusion process must be 
taken into account to properly describe the dynamical evolution of the neutrino asymmetry. 



IV. REGION WITH NO OSCILLATIONS 

In Ref. it was shown that in the framework of the static approximation the sign of 
the final value of total lepton number is the same as the sign of baryon asymmetry term, L, 



and this is due to the presence of a correcting term into the equation (see also Ref. for 
a detailed discussion). Furthermore, if the initial value of the effective total lepton number 
is the same as that of L, then the effective total lepton number never changes sign during 
its evolution [otherwise it changes sign only once]. We will simply refer to this situation as 
no oscillations. 

In Ref. it was also shown that for small mixing angles and |(5m^| ~ 10~^eV^, one 
expects the static approximation to be valid. In this section we present the results of a 
systematic research in the space of mixing parameters that confirm this expectation. This 
research has been done solving numerically the QKE. The repopulation term has been taken 
into account in two different ways. 

In the first one we simply describe the active neutrino distribution assuming a thermal 
distribution (/^^ = /eg). This assumption is equivalent to saying that the collisions are 
able to instantaneously refill the states depleted by the oscillations. Such an assumption is 
justified considering that for temperatures T 3> IMeV (which is typically the case since the 
region where the sign is determined is during the period of exponential growth which occurs 
for Te ~ 3 MeV for 6m'^ ~ 10"^ eV^) and values of < 10-^ one has respectively T :$> H 
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and r ^ ^ua^^^- In this case in the equations (H) it is approximately vahd to replace Pq 
and Pz with the neutrino distributions Za = feqifo and Zg = fus/ foy via: 

Pq = Za + Zg, Pz = Za — Zg. (21) 

The differential equation for P^ gives a differential equation for Zg 

^ = -^f^(^)Py(^)- (22) 

The differential equation for Pq becomes redundant because simply: 

dPo _ dzs dza , . 

with Za = feq/fo, fully determined by the thermal equilibrium assumption |. 

In the second way we calculated the repopulation term using the repopulation equation 
in Eq.(P). The comparison shows nice agreement, as illustrated in fig. 4a,b. Figure 
4a is for the parameter choice dm? = —1 eV^, sin^ 26'o = 10~^, while figure 4b is for 
Sm"^ = —10 eV^, sin^ 29q = 10"''. The solid line is the numerical integration of the QKEs, 
Eq-® while the dashed line assumes the thermal active neutrino distribution (discussed 
above). Note that the two computations were done quite independently using different 
codes. Also, two different initial values of L^'^^ were chosen along with a slightly different 
initial temperature. We have also done many other examples which we do not show for 
environmental reasons (i.e. to save trees). This confirms the expectations about the validity 
of a thermal equilibrium assumption. Note that we only show the high temperature evolution 
since this is the region where the sign is determined. (Recall that the final magnitude of the 
neutrino asymmetry is in the range 0.23 ~ L^^ ~ L'^"-^ /2 ~ 0.37 as illustrated in figures 
la, 2 a and is discussed in detail in Ref. 0). 

For fixed values of |5m^|, we studied the evolution of lepton number for increasing values 
of the mixing angle. In the figure 5 we show an example = —10 eV^. It appears evident 
that for the shown values of mixing angles, the sign of lepton number does not oscillate. As 
we will discuss later on, for small mixing angles this result is expected. Of course changing 
the initial conditions does not change this result as we illustrate in figure 6 Q. This figure 



^ A naive interpretation of the equation (P) for Pq, would suggest that in the limit of thermal 
equilibrium dPQ/dt = 0. This result is clearly incorrect because on the contrary it corresponds 
to the case that interactions are off and repopulation is not active at all. The assumption of 
thermal equilibrium physically corresponds to the case that the rate T{x) of the collisions, that are 
refilling the states of momentum x, is so strong that the difference f^^ — fy^^ is very small. Thus 
it corresponds to the limit F ^ oo with dPo/dt finite given by the Eq. (|23|), from which it follows 
that 



^ Of course qualitatively different behaviour does occur if the initial value of L^^ is large enough, 
roughly, i^^"***"' ~ 10^^ [21,23|. For such a large initial value, the oscillations cannot effectively 



destroy L^"^) and L^"^) remains of order L™***"' until lower temperatures where it eventually further 



increases. 
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considers Vr ^ oscillations with Sm'^ = —10 eV"^ and sin^ 2^0 = 10~^ with four different 
initial values of L^^ employed. 

The result of our careful study of the sign of the neutrino asymmetry is given in fig. 7. 
As the figure shows, the parameter space breaks up into two regions, a no oscillation region 
and an oscillation region. Looking at the no oscillation region, one can easily identify two 
different parts, one part at low sin^ 29q: 

sin^ ~ 10"^ (24) 

and one part at large sin^2^0; 

sin^ 260 (5m2 ~ 3 X 10"^ eV"^. (25) 

The existence of the no oscillation region at low sin^ 29q is consistent with the physical 
insight provided by the static approximation approach. Recall that in this approximation, 
oscillations of sign do not occur . Furthermore this approximation is valid for sin^ 29q 
sufficiently small. This is because the rate of change of lepton number during the exponential 
growth epoch increases proportionally with the value of sin^ 26q. 

In the region of large 3 x 10~^eV^^/(5m^ ~ sin^ 26'o ~ 10~^ oscillations evidently arise at 
the critical temperature. This is not inconsistent with the static approximation, since it is 
not valid in this region. At the moment we do not have a clear analytic description of the 
onset of rapid oscillations, since the QKE's themselves do not provide much physical insight. 
We can give however a qualitative partial explanation noticing that for a <^ b [recall that 
a, b are defined in Eq.(|^)], 

~50sin2^o. (26) 

Thus, the line sin^ 26*0 ~ constant (10~^) also corresponds to P/D constant (10~^) at 
the resonance. At the resonance the quantity (3/D gives approximately the value of the 
ratio between the mean interaction length and the matter oscillation length. Until this 
quantity is much less than unity the coherent behaviour of neutrino oscillations at the 
resonance is averaged out by the effect of collisions and the static approximation provides 
a good description for the evolution of the neutrino asymmetry. When this quantity starts 
to approach unity this starts to be not true any more. Thus, during the initial exponential 
growth, it is possible that some small effect due to the oscillations between collisions, which 
is an effect not included in the static approximation, may be responsible for the rapid 
oscillations of the neutrino asymmetry. 

Of course one may wonder why the oscillation region stops when the mixing angle be- 
comes large enough, Eq. (pSj) . Actually there is quite a simple explanation for this result. 
For large mixing angles such as these, it is known that the exponential generation of lepton 
number is delayed because of the sterile neutrino production At the same time, the 

initial exponential growth is considerably weaker. A consequence of the initial slow expo- 
nential growth of lepton number the static approximation is valid and oscillations in the sign 
of lepton number cannot occur while the neutrino asymmetry is growing relatively slowly. 
Thus this seems to explain why the no oscillation region can extend also to large mixing an- 
gles. Infact, in fig. 8 we show how for \5m?\ = 10^-^ eV^ the critical temperature is delayed 
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and the rate of growth decreased just when the oscillations would be expected to appear at 

~ 10~^. This is surely the reason why for |5m^| ~ 10^'^ eV^, the oscillations do not occur 
for any value of mixing angle. On the other hand for a fixed value of |(5m^| ~ 10^' ''eV^ the 
oscillations can occur between ~ 10~^ and ~ 3 x 10""^ (eV^/|5m^|). 

We have now to discuss whether the presence of statistical fluctuations with the ampli- 
tude estimated in section |T|, can give any effect in the "no oscillations" region. 

In the static approximation framework, one can show the existence of an approximate 
fixed point, that for T > Tc is stable and drives the effective total lepton number - L*^"^ to 
values of order 10"^^ — 10~^^ prior to the onset of the instability (see [0). This is clearly 
confirmed looking at the solutions in figures 1-6 and 8, obtained solving directly the QKEs. 

Therefore, from our estimation of the maximum amplitude of the statistical fluctuations, 
Eq. (^O]) , we can safely exclude any influence of them in determining the final sign of lepton 
number. We can thus conclude that for values of mixing parameters in the region where the 
lepton number does not oscillate, given a certain fixed value of L, the final sign of lepton 
number is the same as that of L and cannot randomly change in different spatial points of 
the early universe. 

In fig. 7 we have also showed the contour (dot-dashed line plus the dotted line) of 
the allowed region for the solution z/^ — Ug to the atmospheric neutrino anomaly, in the 3 
neutrino mixing mechanism proposed in and reanalyzed in 0,^]. Even assuming that an 
overproduction of sterile neutrinos would occur in the parameter region of rapid oscillations, 
as claimed in it appears evident that the mechanism can still allow acceptable values of 
the tau neutrino mass and is not inconsistent with standard BBN. 

We want now to investigate more in detail the nature of the rapid oscillations. 



V. REGION OF RAPID OSCILLATIONS 

The study of the rapid oscillation region is much harder from a numerical point of view. 
It means that to describe it properly the numerical account of momentum dependence must 
be much more accurate (on this point more details are given in the numerical Appendix). 
In particular a much higher resolution in momentum space is necessary to describe the 
distributions. This makes difficult a systematic exploration of this region for different values 
of mixing parameters. 

We mainly performed the runs for 5m? = — lOeV^, a particularly interesting value for 
the three neutrino mixing mechanism mentioned at the end of previous section. In fact it 
would correspond to minimum values of the tau neutrino mass included in the allowed region 
and compatible with structure formation arguments. 

The transition from the region with no oscillations to that one with rapid oscillations is 
very sharp, meaning that the width of the boundary between the two regions is very narrow. 
In figure 9 we give an example of this transition for a fixed value of Srn^ = — lOeV^. 
Changing the value of log(sin^^o) for a very small quantity, ~ 0.1, is sufficient to pass from 
one region to another. 

The oscillations take place in a limited interval of temperatures. This interval increases 
going from points at the border to points in the middle of the rapid oscillations region in the 
mixing parameter space. This feature again confirms the expectation that the oscillations 
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are the effect of the breaking of the adiabaticity condition during some period of the lepton 
number growth and that moving toward the center of the rapid oscillations region, this 
period becomes longer. 

A first observation about the rapid oscillations is that they develop with a temperature 
period of about (10^^ — 10^^) MeV, depending on the value of |5m^|. This period is roughly 
described by the interval of temperature corresponding to the interaction mean time of 
neutrinos F"^ for x ~ 2.2: 

A„,r ^ 4? . 2 MeV (Mp^)^ 5 X 10- {^^f MeV. (27) 

This is in agreement with the expectation for which the coherent nature of the oscillations 
is damped by the collisions unless the effective potential is rapidly changing. 

In the numerical appendix we will show how increasing the resolution in the momentum 
space for the description of the distributions, has the effect of diminishing the amplitude of 
the oscillations. For the example in figure 9 with sin = 10-6•^^ this effect is such that 
if the resolution is not high enough, then the amplitude of the oscillations is so large that 
the total lepton number changes sign, while when a good resolution is employed and the 
numerical solution becomes stable with changing the resolution, the sign changes disappear. 
If the oscillations are artificially amplified by a numerical error, this can induce the suspect 
that the sign changes are not a real feature of the solutions. However for a choice of mixing 
parameters more inside the rapid oscillations region, like already for sin^ = 10~^'^ in 
figure 9, it seems that the sign changes do not disappear increasing the accuracy. As we 
will say in more detail in the appendix, on this point we cannot be conclusive because we 
cannot reach the required resolution to have perfect numerical stability. It seems that more 
advanced computational means are necessary to get a definitive conclusion, while in this 
work we prefer to say that there is an indication that the lepton number oscillations change 
the sign of the solution during the growth regime. We think however that on this point it 
would be desirable to have some rigorous analytical demonstration, that in our opinion is 
still missing, even though some attempts have been recently done WA . 



Another important problem is whether the final sign of lepton number is sensitive or not 
to the small statistical fluctuations necessarily present in the early Universe, as we saw in 
the third section. In our numerical calculations we observe that changes in the initial lepton 
number are able to change the final sign with a choice of mixing parameters well inside 
the rapid oscillations region where a large number of oscillations take place. While moving 
toward the border this effect tends to disappear. This behaviour seems to confirm an effect 
of 'dephasing' of the solutions as described in ||T^. This means that even though solutions 



with different initial conditions converge at the same values prior the onset of the instability 
at the critical temperature, (as shown for example in the figure 6), small relic differences 
will afterwards develop, resulting in growing phase differences during the oscillatory regime 
and in a possible final sign inversion. However also on this point we cannot be conclusive. It 
is infact necessary first to get a perfect numerical stability in the description of the neutrino 
asymmetry oscillations and then to prove that statistical fluctuations as tiny as those present 
in the early Universe are able to yield a dephasing effect. For this second step it is then 
necessary that the numerical error is lower than the statistical fluctuations. Thus we cannot 
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be conclusive because we cannot reach such a computational accuracy. This means also that 
our analysis cannot neither exclude or demonstrate the possibility of a chaotical generation 
of lepton domains, but in any case we can say that this possibility is limited to a region 
of mixing parameters contained in the rapid oscillations region, in a way that the neutrino 
asymmetry undergoes a number of oscillations high enough that two initially close values 
of the neutrino asymmetry, result, in the end, in a different final sign. The possibility 
of drawing the exact contour of such a region, inside the rapid oscillations one, is again 
something that requires a numerical analysis beyond our present computational means. 

We can thus conclude this section saying that the expectation of the appearance of 
coherent effects on the evolution of lepton number when an adiabatic condition breaks is 
confirmed by the numerical analysis. This suggests that the oscillations are not just an 
effect of the numerical error but an intrinsic feature of the equations. The amplitude of 
these oscillations is much likely large enough to change the sign of the neutrino asymmetry, 
but on this point we cannot be conclusive, because our maximum accuracy in the description 
of the momentum dependence is not sufficient to get a good numerical stability and we prefer 
to speak of indication of sign changes. The demonstration that the changes of sign can induce 
a chaotical generation of lepton domains is beyond our numerical analysis, but in anycase 
this can occur only in a region confined within the region of rapid oscillations shown in figure 
7. Therefore any application that makes use of a chaotical generation of lepton domains as 
a general consequence of the QKE, is not justified. 

VI. CONCLUSIONS 

Large neutrino asymmetry is generated in the early Universe by ordinary - sterile neutrino 
oscillations for a large range of parameters. We have made a numerical study of the final 
sign of this asymmetry and our results are summarized in figure 7. In the space of mixing 
parameters we identified a no oscillation and an oscillation region. In the no oscillation 
region we could show perfect agreement with the predictions of the static approximation for 
temperatures T ~ at which the MSW effect is inhibited by the presence of the collisions. 
Moreover perfect numerical stability of the solutions was obtained. These two facts do not 
leave any room for doubt about the generation of a neutrino asymmetry and its final value. 

We have also estimated an upper limit [see Eq.(pOD] for the size of statistical fluctuations 
of the lepton number which is very tiny. In view of this, for a choice of mixing parameters in 
the no oscillation region, the final sign of the neutrino asymmetry is the same in all points 
of space unless there happens to be larger fluctuations (of non-statistical origin) present. 
In this way the possibility of a chaotical generation of lepton domains is excluded in the 
no oscillation region. One consequence of this is that the three neutrino mixing mechanism 
(proposed in [^) allowing the i/^ ^ Ug solution to the atmospheric neutrino anomaly to be 
consistent with a stringent BBN bound of N^J^j^ < 3.5 is still viable. 

We have also been able to plausibly justify the existence of the rapid oscillations region 
in terms of a breaking of the static approximation during the lepton number growth. 

The study of the rapid oscillations in the oscillation region encounters many numerical 
difficulties. The oscillations seem to have an amplitude sufficiently large to change the sign 
of the total lepton number during the growth regime. We however think that further analysis 
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are required to be conclusive on this point. We also think that at the present the numerical 
analysis cannot neither exclude or demonstrate the possibility of a chaotical generation of 
lepton domains. This is because the error introduced by the numerical analysis is much 
larger than the tiny statistical fluctuations present in the early universe. In anycase this 
possibility is confined to a special choice of the mixing parameters, corresponding to some 
region contained within the region of rapid oscillations. 

Finally we should also mention that while we have focussed on i/q, oscillations 
in isolation, our results will be applicable to realistic models with sterile neutrinos. In the 
more general case of three ordinary neutrinos and one or more sterile neutrinos, the neutrino 
asymmetry generation occurs first for the oscillations with the largest |5m^| (with Sm^ < 0). 
The sign of this asymmetry (which we assume to be L,^^ for definiteness) will be the same as 
L except possibly in the oscillation region where it may depend on (5m^, sin^ 29o. The other 
oscillations with smaller |5m^| may also generate significant neutrino asymmetry, although 
this generation occurs latter on (see e.g. Ref. for some examples). The sign of these 
asymmetries will ultimately depend only on the sign of L^,^ since this dominates L for these 
oscillations. 

NUMERICAL APPENDIX 

In this Appendix we want to describe in more detail our numerical procedure for solving 
the QKE's. As we said already in the text we employed two different codes, one using 
the thermal equilibrium assumption [code A), the other using the expression (|^) for the 
repopulation account {code B). The code B has been already described in [Q and we refer 
the reader to this reference for details. Here we describe common features of the two and 
specifications of the code A. 

We first discuss the time step and then the momentum integration is described. 

The time step of integration is adjusted in a way that it is halved until the required 
accuracy (the relative error on all the set of variables of the system) e is reached. By this 
we mean that for each variable we introduce a time step At. At every time step during 
the evolution, say t = we compute Z{tx + At) and compare this with the results of 
integration with a half step size, symbolically, Z[tx + At/2 + At/2). Explicitly, the step size 
is halved until every integration variable (which we have denoted collectively as Z) satisfies, 

\Z{t^ + At) - Z(t^ + At/2 + At/2)| < e\Z{t^ + At) + Z{t^ + At/2 + At/2)|/2. (28) 

The above procedure is quite standard, and is called Merson' s procedure' or step dou- 
bling p4| , p5| and together with a fourth order Runge-Kutta solver is considered the most 
straightforward and safe technique 0. 



° An evolution of this procedure is given by the "Fehlberg" or embedded Runge-Kutta method. 
This method should be twice as fast as Merson's one and although it should be safe enough, we 
have adoped the more conservative Merson's procedure for code A. Another possibility is offered by 
the Bulirsch-Stoer method that is by far the most efficient but it has many cautions for its usage. 



14 



We now discuss the discretization of momentum and the integration around the reso- 
nances procedure. The range of values of momenta (p/T) is [xi,X7v] which is discretized on 
a logarithmically spaced mesh: 



Xi X 



lOO-D' 



J 



1,2, 



N 



(29) 



where N is the number of bins and A = log(xiv/a;i)/(A^ — 1). The initial temperature is 
chosen in a way that the initial resonant momentum is 

The resonance momentum for neutrinos, Xres = Pres/T can be obtained by solving A(x) = 
[with A(x) defined in Eq.(|ll])] and is given by 



2Xi ^ \ \2Xj ^ Xi 



fX2 



(30) 



where c = cos 2^0 ^^^d Xi = b/x'^,X2 = a/x [recall a,b are defined in Eq.(T^)]. Note that 
Xi 2 are independent of x. The resonance momentum for antineutrinos can be obtained by 
replacing X2 — >■ —X2 in the above equation. The width of the resonance is given by the 
following expression f\ 



6x = 2 



sin^ 200 + dl 62 



c + 6 



= 5r 



res ■^resj 



(31) 



where we introduced the logarithmic width of the resonance 6res = Slnx and do = 
2presD/b6m'^. Numerically, do ~ 2(0.8) x 10"^ for a = /i, r (e). Note that 6res is roughly 
independent of T for T ~ Tc (and is of order do) while for T <^ Tc, 5res decreases to ~ sin 26o. 
It is also roughly independent of the momentum x and this is the reason for which it is more 
convenient to use a logarithmically spaced mesh rather than a linearly spaced one. 

At each step of integration we calculate both 
jres with 



:.res and (5res and we discretize them: Xres 



Jr 



Int 



1 7 2^res 

A xi 



+ 1 



(32) 



and we define Aj = 1 + Int((5rcs/A). Considering that Sj-cs = S{lnx) ^ 2.3 x Slogx, this 
means that Aj corresponds roughly to twice the width of the resonance in a logarithmic 
scale and in units of A. At each time step we integrate only on those bins among [1, N] in 
the interval [jmin, jmax] where: 



Same considerations hold for predictor-corrector methods that are considered somehow worse than 
all of the previous ones. As we had to study unexplored features of still 'young' equations, we 
preferred to adopt the step doubling method (see [^ for further details). 



^This is derived in Ref. [^. The form given in Eq.(31) is simpler, but mathematically equivalent 
to A/pres where A is given in Eq.(28) of Ref. Q. 
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Jmin Jrcs P ) Jmax Jrcs ~l~ P Aj 



(33) 



In this way we are integrating symmetrically around the resonance (we are taking equal 
numbers of bins with lower and higher momenta), but logarithmically, while in the code B 
this is done linearly [Q. Clearly jmin and jmax are constrained within [1, A^]. This procedure 
is done at each step so that the interval of integration follows dynamically the resonance. 
Initially the neutrino and anti-neutrino resonances approximately coincide due to the ap- 
proximate fixed point L^") ^ 0. However during the exponential growth the resonance splits 
in two, one at low momenta for the antineutrinos and one at high momenta for the neutrinos 
if L(") > or vice versa if < 0. The split occurs gradually: first the initial interval 
enlarges and at certain point it splits, so that eventually the two intervals do not overlap. 

Let us introduce the number of bins which we integrate over, which for neutrinos we 
denote by and is given by = min(A^, Jmax) — max(l, Jmin) (and similarly for antineu- 
trinos we introduce the quantity Mp). So at each step of integration we solve a system of 
3(M^ + Mp) + 1 equations in the code A and A{M^ + Mp) + 1 in the code B. 

In this way one has 6 numerical parameters: e,xi,X]s[,Xin, p, N. One has to check care- 
fully that the numerical solution is 'stable' changing these parameters in the direction of 
greater accuracy. 

For the parameter e we checked that e = 10~^ is a safe value. A safe choice for (xj„, Xi, x^) 
is (0.1, 10"'^, 20) [xin = 0.1 corresponds approximately to the choice Tj„ ~ 3 x Tc where Tc 
is given by Eq. ([T8|)]. 

The minimum value for A^ to get an accurate description, depends on the values of the 
mixing parameters, whether they are outside or inside the rapid oscillations region. The 
number of bins A^, necessary to have a stable solution, is that one for which the number of 
bins within a width of resonance, Aj, is large enough. A value Aj ~ 10, that corresponds 
to A^ ~ 2^^, is sufficient in the no oscillatory region to have stable solutions as shown in 
figure 10a. 

A first observation concerning the value of p to be chosen is that p ~ 100 corresponds 
roughly to integrate on all the range of momenta 10^'^ — 20 0. For values of mixing 
parameters in the no oscillation region, p = 4 is a good choice in order to study the problem 
of sign. In figures 10a, 10b we present the result of a check that shows the numerical 
stability of our results for the indicated values of the mixing parameters. In figure 10a 
the mixing parameters correspond to a point close to the border with the rapid oscillations 
region. 

More precisely the minimum value of p to get the stability, in the no oscillations region, 
depends on the initial value of the total lepton number. If one choose L-"^ = 0, a relative 
low value of p = 4 is sufficient to get a full stability at all temperatures in an impressive way. 



we start from Xj„ = 0.1 and if we consider the case a = /i,r so that 6res — do — 0.02 (for 
T ~ Tc), then the choice p = 100 is equivalent to integrating on the interval 10^'^ < x < 10. At 
the critical temperature Xres ~ 2, the integration interval moves at values 0.02 < x < 20. This 
means that the integration around the resonance is not necessarily more approximated than the 
case when a static interval of integration is chosen, it is simply a way to select a reduced interval 
of momenta but dynamically, taking into account the resonant behaviour of the process. 
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On the other hand if one starts from non vanishing values of the initial lepton number, the 
approximation of integrating around the resonance gives an error at high temperatures, in 
the regime during which the lepton number is destroyed. This is also clearly shown in figure 
10a and 10b for the indicated mixing parameters. This fact has already been described and 
justified in the main text. 

In the oscillatory region, in order to get stability, the necessary number of bins inside a 
resonance width seems to hugely increase upto about Aj ~ 500, obtained for ~ 2^^, as 
shown in figures 11a, lib. It is remarkable to notice how an inaccurate description of the 
momentum dependence amplifies the amplitude of the oscillations and for sin^ 29q = 10~^'^^ 
it makes the solution changing sign. However for an higher value of the mixing angle the 
sign changes are present also when the maximum accuracy is employed. The stability of 
the solutions in the oscillatory regime is clearly not as good as outside this regime, even 
though there is still some indication on the behaviour of the solutions. Not only it would be 
desirable to make checks increasing the number of bins N, but also increasing the width of 
the region of integration around the resonance, that means the parameter p. Unfortunately 
this is not possible not only because the required computational time becomes too long, but 
also because the dimensions of the arrays employed to describe the distributions requires 
an amount of RAM memory than our machines do not have. It would be also desirable to 
make checks increasing the numerical precision, in our case passing from a double precision 
(the round off error is about 10^^^) to a quadruple precision. This because at the onset of 
the instability the rate of growth of neutrino asymmetry is the difference of two opposite 
quantities that are much higher. Unfortunately we do not have at the moment the possibility 
to perform quadruple precision runs, also because again the time for the runs becomes too 
long. For all these reasons we prefer to speak of a strong 'indication' that the oscillations 
are able to change the sign of the neutrino asymmetry. 

In the end of this appendix we want to mention that also for an accurate evolution at 
low temperatures T <^ Tc, large values of are required (depending on sin2^o)- This is 
because, the resonance width becomes quite narrow at low temperatures [see the discussion 
above, following Eq.(^)]. 

This appendix has been stimulated by the scepticism expressed in [|I0| about the numer- 
ical procedure employed in the previous existing works dealing with the numerical solutions 
of the exact QKE. It should be clear now that the results and the procedure were already 
rigorous at that time, just there was not a compelling necessity to provide all the tedious 
details as we do now. 

Acknowledgements 

We thank Nicole Bell, Roland Crocker, Ray Volkas and Yvonne Wong for valuable dis- 
cussions. P. Di Bari wishes to thank Maurizio Lusignoli for helpful discussions and for his 
support and encouragement to the collaboration with R. Foot and the Melbourne University. 
He also wants to thank Piero Cipriani, Maria Di Bari, Maurizio Goretti, Trevor Hales and 
Rob Scholten for help with his portable computer making his work in Melbourne much easier 
and faster. R.F. acknowledges the kind hospitality of Paolo Lipari and Maurizio Lusignoli 
at the University of Rome during a visit where this work was initiated. 



17 



REFERENCES 



[1] R. Foot, M. J. Thomson and R. R. Volkas, Phys. Rev. D53, 5349 (1996). 

[2] R. Foot and R. R. Volkas, Phys. Rev. D55, 5147 (1997). 

[3] R. Foot and R. R. Volkas, Phys. Rev. D56, 6653 (1997); Erratum-ibid, D59, 029901 
(1999). 

[4] R. Foot, Astropart. Phys. 10, 253 (1999). 

[5] N. F. Bell, R. R. Volkas and Y. Y.Y.Wong, Phys. Rev. D59, 113001 (1999). 

[6] P. DiBari, P. Lipari and M. Lusignoh, |hep-ph/9907548 . 



[7] For BBN implications of the neutrino asymmetry in various 4 and 6 neutrino models, 
see Ref. and also N. F. Bell, R. Foot and R. R. Volkas, Phys. Rev. D58, 105010 

(1998); R. Foot, |hep-ph/9906311| (to appear in Phys. Rev. D). For applications to models 
with mirror neutrinos, see R. Foot and R. R. Volkas, Astropart. Phys. 7, 283 (1997); 
|hep-ph/9904336| (to appear in Phys. Rev. D). 

K.Enqvist, K.Kainulainen and J.Maalampi, Nud.Phys. B349 754 (1991); R.Barbieri and 
A.D. Dolgov, Nud.Phys. B349 743 (1991). 



[9 
[10 

[11 

[12 

[13 
[14 
[15 

[16 
[17 
[18 
[19 
[20 
[21 
[22 
[23 
[24 
[25 



D.P. Kirilova and M.V. Chizhov, Nud.Phys. B 534, 447 (1998). 
A.D. Dolgov, S.H.Hansen, S. Pastor and D.V.Semikoz, [hep-ph/ 99 10444 . 
X. Shi and G. Fuller, Phys. Rev. Lett. 83, 3120 (1999). 
X. Shi, Phys. Rev. D54, 2753 (1996). 

K. Enqvist, K. Kainulainen and A. Sorri, Phys. Lett. B464, 199 (1999). 

A. Sorri, |hep-ph/99H"366| . 

R. A. Harris and L. Stodolsky, Phys. Lett. 116B, 464 (1982); Phys. Lett. B78, 313 
(1978); A. Dolgov, Sov. J. Nucl. Phys. 33, 700 (1981). 

B. H. J. McKellar and M. J. Thomson, Phys. Rev. D49, 2710 (1994). 

L. Stodolsky, Phys. Rev. D36, 2273 (1987); M. Thomson, Phys. Rev. A45, 2243 (1991). 

K. Enqvist, K. Kainulainen and M. Thomson, Nucl. Phys. B 373, 498 (1992). 

L. Wolfenstein, Phys. Rev. D17, 2369 (1978). 

D. Notzold and G. Raffelt, Nucl. Phys. B307, 924 (1988). 

K. Enqvist, K. Kainulainen and J. Maalampi, Phys. Lett. B244, 186 (1990). 

P. Di Bari, |hep-ph/99112Ti[ 

R. Foot and R. R. Volkas, Phys. Rev. Lett. 75, 4350 (1995). 

G. N. Lance, 'Numerical methods for high speed computers' (London 1960). 

W.H. Press, S.A.Teukolsky, W.T.Vetterling and B.P. Flannery, 'Numerical Recipes, the 

Art of Scientific Computing', Cambridge University Press (1994). 



FIGURE CAPTIONS 

Fig. 1 Evolution of the effective total lepton number, L'^'^\ for Ur — Us oscillations for the 
parameter choice, 6m'^ = —10 eV"^, sin^ 29q = 10"''. In Fig. la, the solid line is the numerical 
solution of the QKEs, Eq.(^ integrating around the resonances using / = 12 resonance 
widths and in a linearly symmetric way (the 'code B' has been used: see the Appendix and 
Q for its definition), while the dashed line is the static approximation (described in the text) 
also integrating around the resonance. Fig. lb is a comparison of the static approximation 
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integrating around the resonance (solid line) with integrating over the entire momentum 
range 0.01 < p/T < 20 (dashed line). Note that L = 5 x 10~^° is used and the initial value 
of is zero in this example. 

Fig.2a,b Same as figures la,b except that the oscillation parameters are changed to 
= -100 eV^, sin^2^o = 10"^ 

Fig. 3 Numerical solution of the QKE's (solid line) integrated over the entire momentum 
range 0.1 < x < 12 for the parameter choice = —0.01 eV^, sin^ 29o = 10~^. The initial 
value of is zero in this example. Also shown is the static approximation integrated over 
the same momentum range and same initial conditions (dashed line). 

Fig. 4a,b High temperature evolution of the effective total lepton number, L^'^\ for 
Ur — Ug oscillations. Figure 4a is for the parameter choice, dm? = —1 eV^, sin^ 2^o = lO^*'^ 
(with L^l^itiai = L = 5 X 10"^°, i.e. initial = 0) while figure 4b is for Sm^ = —10 eV"^, 
sin^ 26^0 = 10'^ (with L^lniUai = 10~^^). The solid line is the numerical solution of the QKEs, 
Eq.(^, while the dashed line is the solution of the QKEs assuming a thermal distribution 
for the active neutrino (as discussed in the text). 

Fig. 5 An example of the high temperature evolution of the effective total lepton 
number, L^^) for dni^ = -10 eV^ for sin^ 2^o = 10"^ (dotted line), 10"*^ (dashed line) 10"^ 
(solid line). 

Fig. 6 An example of the high temperature evolution of the effective total lepton 
number, L^'^^ for Sm^ = —10 eV'^ for sin^ 29q = 10^^ with the initial values of L^^ given 
by L^^ = 10^^° (solid line), 10~^ (long dashed line), 10~^ (short dashed line), 10~^ (dotted 
line) . 

Fig. 7 Region of parameter space where the final sign does and does not oscillate 
for z/t- ^ oscillations. Also shown is the 'allowed region' (which we define below) for the 

^ maximal oscillaton solution to the atmospheric neutrino anomaly. This is the region 
where the L,^^ is generated rapidly enough so that the sterile neutrino is not equilibrated 
by either ^ Vg maximal oscillations (region above the dashed-dotted line) or by v^- ^ Vg 
oscillations (region to the left of the dotted line). 

Fig. 8 The evolution of the total lepton number for a choice of mixing parameters just 
outside the top of the rapid oscillations region. The sterile neutrino production increases 
with the mixing angle. It appears clear the double effect of the sterile neutrino production 
in delaying the onset of the growth, lowering the critical temperature and in depressing the 
rate of growth of lepton number. The latter prevents the arise of rapid oscillations when the 
mixing angle increases. 

Fig. 9. Transition from the no oscillations region to the rapid oscillation region increas- 
ing the value of the mixing angle for a fixed value of 5m? = — 10 eV^. 

Fig. 10 Check of numerical stability for two different choices of mixing parameters inside 
the no oscillation region for 5m? = —10 eV^ and sin^ 26o = 10~^-^ (a), sin^ 26*0 = 10"'' (b). 
It appears clear the validity of the resonant approximation in this region when = 0. On 
the other hand for L-"^ = L = 5 x 10~^^, the resonant approximation produces an error in 
the regime of high temperatures (T ^ T^) when the initial lepton number is destroyed. 

Fig. 11 Check of numerical stability in the oscillatory regime for = — 10 eV^ and 
sin^ 200 = 10-^-^^ (a), sin^ 2^o = lO"^'^ (b). 
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figure 8 
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figure 9 
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figure 10a 
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figure 11a 
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figure 1 1b 
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